# Code for Table M7: Clustered standard errors, magistrate × year × month × case severity

rm(list = ls())

## ---------------------------------------
## Load Packages 
## ---------------------------------------

library('AER')
library('ivpack')
library('dplyr')
library('data.table')


## ---------------------------------------
## Load Data and Functions
## ---------------------------------------
load("../data0812.RData")

directory <- "../functions/"
functions <- list.files(directory)  
loadfunctions <- sapply(functions, FUN = function(x)source(paste0(directory, x)))

## ---------------------------------------
# Construct Instrument
## ---------------------------------------
data0812 <- constructIV(data0812)

## ---------------------------------------
#  Select Last Case Before Election
## ---------------------------------------
data0812 <- lastCase(data0812)

## ---------------------------------------
## Two Stage Least Squares Specifications
## ---------------------------------------
time.controls <- "as.factor(court_time1) + as.factor(court_time2) + as.factor(court_dow) + as.factor(court_shift) + as.factor(totOGS2)"
case.controls <-   "as.factor(any_drug_2) +  as.factor(any_violent_2) + as.factor(fire_arms_2) +  as.factor(any_rob_2) + as.factor(any_dui_2) + as.factor(prior_offender_2)"
demo.controls <- "age_2012 + I(age_2012^2) + Female + as.factor(race) + vote2008 + as.factor(noteli08) + regis_before"

outc.1 <- "vote2012"
endo.1 <- "pti"
inst.1 <- "judgeiv"

form.1 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "|", inst.1, "+", time.controls))
form.2 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "|", inst.1, "+", time.controls, "+", demo.controls))
form.3 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "+", case.controls, "|", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls))

m1a1 <- ivreg(form.1, data = data0812)
m1a2 <- ivreg(form.2, data = data0812)
m1a3 <- ivreg(form.3, data = data0812)


data0812[, clu_id := .GRP, by = c("judge_cat", "court_time1", "court_time2", "totOGS2")]

m1a1_c <- round(coeftest.cluster(data0812, m1a1, cluster1="clu_id")["pti", 1:2],3) 
m1a2_c <- round(coeftest.cluster(data0812, m1a2, cluster1="clu_id")["pti",1:2 ],3)
m1a3_c <- round(coeftest.cluster(data0812, m1a3, cluster1="clu_id")["pti",1:2 ],3) 
row1 <- cbind(m1a1_c, m1a2_c,m1a3_c)
row1n <- rep(m1a3$nobs,3)

## ---------------------------------------
## Row 2: Racially Disparate Impact
## ---------------------------------------

outc.1 <- "vote2012"
endo.1 <- "pti * as.factor(race)"
inst.1 <- "judgeiv * as.factor(race)"

form.1 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "|", inst.1, "+", time.controls))
form.2 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "|", inst.1, "+", time.controls, "+", demo.controls))
form.3 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "+", case.controls, "|", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls))

data0812$race <- as.factor(data0812$race)
data0812 <- within(data0812, race <- relevel(race, ref = "Black"))

data0812[, clu_id := .GRP, by = c("judge_cat", "court_time1", "court_time2", "totOGS2")]

m1a1 <- ivreg(form.1, data = data0812)
m1a2 <- ivreg(form.2, data = data0812)
m1a3 <- ivreg(form.3, data = data0812)

m1a1_c <- round(coeftest.cluster(data0812, m1a1, cluster1="clu_id")[c("pti"), 1:2],3)
m1a2_c <- round(coeftest.cluster(data0812, m1a2, cluster1="clu_id")[c("pti"), 1:2],3)
m1a3_c <- round(coeftest.cluster(data0812, m1a3, cluster1="clu_id")[c("pti"), 1:2],3)
row2a <- cbind(m1a1_c, m1a2_c, m1a3_c)

m1a1_c <- round(coeftest.cluster(data0812, m1a1, cluster1="clu_id")[c( "pti:as.factor(race)White"), 1:2],3)
m1a2_c <- round(coeftest.cluster(data0812, m1a2, cluster1="clu_id")[c("pti:as.factor(race)White"), 1:2],3)
m1a3_c <- round(coeftest.cluster(data0812, m1a3, cluster1="clu_id")[c( "pti:as.factor(race)White"), 1:2],3)
row2b <-cbind(m1a1_c, m1a2_c, m1a3_c)

row2<- rbind(row2a, row2b)
row2n <- rep(m1a3$nobs,3)


## ---------------------------------------
## Row 3: Resource Deprivation
## ---------------------------------------
outc.1 <- "vote2012"
endo.1 <- "pti * as.factor(quants_inc)"
inst.1 <- "judgeiv * as.factor(quants_inc)"


form.1 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "|", inst.1, "+", time.controls))
form.2 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "|", inst.1, "+", time.controls, "+", demo.controls))
form.3 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "+", case.controls, "|", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls))

data0812[, clu_id := .GRP, by = c("judge_cat", "court_time1", "court_time2", "totOGS2")]
data0812 <- data0812[!is.na(quants_inc), ]

m1a1 <- ivreg(form.1, data = data0812)
m1a2 <- ivreg(form.2, data = data0812)
m1a3 <- ivreg(form.3, data = data0812)
m1a1_c <- round(coeftest.cluster(data0812, m1a1, cluster1="clu_id")[c("pti"),1:2 ],3)  
m1a2_c <- round(coeftest.cluster(data0812, m1a2, cluster1="clu_id")[c("pti"), 1:2],3)
m1a3_c <- round(coeftest.cluster(data0812, m1a3, cluster1="clu_id")[c("pti"),1:2 ],3) 
row3a <- cbind(m1a1_c, m1a2_c, m1a3_c)

m1a1_c <- round(coeftest.cluster(data0812, m1a1, cluster1="clu_id")[c("pti:as.factor(quants_inc)2"),1:2 ],3)  
m1a2_c <- round(coeftest.cluster(data0812, m1a2, cluster1="clu_id")[c("pti:as.factor(quants_inc)2"), 1:2],3)
m1a3_c <- round(coeftest.cluster(data0812, m1a3, cluster1="clu_id")[c("pti:as.factor(quants_inc)2"),1:2 ],3) 
row3b <- cbind(m1a1_c, m1a2_c, m1a3_c)

m1a1_c <- round(coeftest.cluster(data0812, m1a1, cluster1="clu_id")[c("pti:as.factor(quants_inc)3"),1:2 ],3)  
m1a2_c <- round(coeftest.cluster(data0812, m1a2, cluster1="clu_id")[c("pti:as.factor(quants_inc)3"), 1:2],3)
m1a3_c <- round(coeftest.cluster(data0812, m1a3, cluster1="clu_id")[c("pti:as.factor(quants_inc)3"),1:2 ],3) 
row3c <- cbind(m1a1_c, m1a2_c, m1a3_c)

row3 <- rbind(row3a, row3b, row3c)
row3n <- rep(m1a3$nobs,3)


## ---------------------------------------
## Row 4: Resource Deprivation (different baseline)
## ---------------------------------------
data0812$quants_inc <- as.factor(data0812$quants_inc)
data0812 <- within(data0812, quants_inc <- relevel(quants_inc, ref = 2))

m1a1 <- ivreg(form.1, data = data0812)
m1a2 <- ivreg(form.2, data = data0812)
m1a3 <- ivreg(form.3, data = data0812)

m1a1_c <- round(coeftest.cluster(data0812, m1a1, cluster1="clu_id")[c("pti", "pti:as.factor(quants_inc)3"),1:2 ],3)  
m1a2_c <- round(coeftest.cluster(data0812, m1a2, cluster1="clu_id")[c("pti", "pti:as.factor(quants_inc)3"),1:2 ],3)  
m1a3_c <- round(coeftest.cluster(data0812, m1a3, cluster1="clu_id")[c("pti", "pti:as.factor(quants_inc)3"),1:2 ],3)  

m1a1_c <- round(coeftest.cluster(data0812, m1a1, cluster1="clu_id")[c("pti"),1:2 ],3)  
m1a2_c <- round(coeftest.cluster(data0812, m1a2, cluster1="clu_id")[c("pti"),1:2 ],3)  
m1a3_c <- round(coeftest.cluster(data0812, m1a3, cluster1="clu_id")[c("pti"),1:2 ],3)  
row4a <- cbind(m1a1_c, m1a2_c, m1a3_c)

m1a1_c <- round(coeftest.cluster(data0812, m1a1, cluster1="clu_id")[c( "pti:as.factor(quants_inc)3"),1:2 ],3)  
m1a2_c <- round(coeftest.cluster(data0812, m1a2, cluster1="clu_id")[c(  "pti:as.factor(quants_inc)3"),1:2 ],3)  
m1a3_c <- round(coeftest.cluster(data0812, m1a3, cluster1="clu_id")[c("pti:as.factor(quants_inc)3"),1:2 ],3)  

row4b <- cbind(m1a1_c, m1a2_c, m1a3_c)

row4 <- rbind(row4a, row4b)
row4n <- m1a3$nobs


## ---------------------------------------
# Table M7
## ---------------------------------------
cont <- as.data.frame(rbind(row1, row1n,
                            row2, row2n,
                            row3, row3n,
                            row4, row4n))
colnames(cont) <- c("(1)", "(2)", "(3)")
cont$names <-  c("Main result: PTI",
                 "Std. Error",
                 "Obs.",
                 "PTI (Black as baseline)",
                 "Std. Error",
                 "PTI × White",
                 "Std. Error",
                 "Obs.",
                 "PTI (Bottom tercile as baseline)",
                 "Std. Error",
                 "PTI × Middle tercile",
                 "Std. Error",
                 "PTI × Top tercile",
                 "Std. Error",
                 "Obs.",
                 "PTI (Middle tercile as baseline)",
                 "Std. Error",
                 "PTI × Top tercile",
                 "Std. Error",
                 "Obs."
)
cont <- cont[, c("names", "(1)", "(2)", "(3)")]
rownames(cont) <- NULL
cat("\nPrinting Table M7: Clustered Std Errors magistrate-year-month-severity...\n")
print(cont)



